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ABSTRACT 


The consequences of the earth's elasticity are examined for close-earth satellites. 
The ideas of polar motion and earth tides are developed in a form applicable to satel- 
lite studies, since the polar motion, the body tide, and the ocean tide are all suitable 
for study by use of satellites. Analysis of available polar-motion data is performed. 
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INVESTIGATIONS OF EARTH DYNAMICS FROM SATELLITE OBSERVATIONS 

Final Report 


1. INTRODUCTION ; 

The principal aim of this grant, initiated m 1969, was to explore the use of 
precision satellite-tracking data for studying the rotating elastic earth, i. e. , polar 
motion, earth tides, and related phenomena. The data available were camera observa- 
tions taken for geodetic purposes, with a limited amount of precision laser tracking 
data. These tracking data had been used in the determination of the 1969 Smithsonian 
Standard Earth (II) (Gaposchkin and Lambeck, 1971). 

This study had the following aspects: 1) Theoretical investigation of polar-motion 
dynamics; 2) assembly of some mathematical tools; 3) study of the consequences of 
polar motion and tidal deformation on satellite observations and orbits; 4) incorporation 
of these effects into the Smithsonian Astrophysical Observatory (SAO) geodesy programs 
where appropriate; 5) analysis of polar-motion measurements made by classical astro- 
nomical methods; 6) assembly of the available tracking data for polar-motion deter- 
mination; and 7) determination of polar motion based on satellite observations. 

The principal results are as follows: 

1) Improved values of polar motion and tidal parameters can be obtained from 
analysis of satellite-tracking data. Polar motion is most easily determined from 
kinematic effects, and body tides are most easily determined from analysis of satellite 
perturbations. 

2) Analysis of existing tracking data gives reasonable results for tidal parameters, 
but the determination of polar motion was not possible with the data and orbits avail- 
able at that time. 
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3) Polar motion and the motion of the principal axis of the earth can eventually be 
studied through satellite determination of and S^. The orbital effects are approx- 
imately 1 m. 

4) The body tide is most easily studied through the gravitational perturbation of 
satellite orbits. An important contribution to the tidal perturbation of satellite orbits 
comes from the ocean tide. There is so little agreement as to the size of the second- 
degree global ocean tide that it must be considered unknown at this time. Consequently, 
certain components of the ocean tides may have to be determined from satellite obser- 
vations. 

5) The ocean-tide effect has been ignored in previous work on the determination 

of the body tide. Its size could explain the variety of numerical values for the Love 

number k_. 

z 
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2. POLAR MOTION AND EARTH TIDES 


Polar motion and earth tides depend on the elasticity of the earth — more precisely, 

/v> 

on the rigidity p. The theory for polar motion, i. e., free nutation, for a rigid body 
was worked out by Euler. Using the precession to obtain (C-A)/C, Euler predicted a 
period of free nutation of 307 days. Attempts to observe this free nutation were 
unsuccessful until 1890, when Kustner demonstrated the existence of periodic variations 
of latitude. Shortly thereafter, Chandler showed the period for free nutation to be 
approximately 430 days rather than 307. Newcomb then demonstrated that the elasticity 
of the earth would lengthen the period of rigid-body free nutation and estimated the 
rigidity to be approximately that of steel or glass. Poincare found the effects of the 
liquid core would shorten the free period, and Lord Kelvin demonstrated that the oceans 
would lengthen the period. The principal reference on polar motion and its geophysical 
consequences is Munk and MacDonald (1960). 

The theory of forced nutation, assuming a rigid earth, was completed by Oppolzer 
in 1882. Schweydar concluded that the elasticity would not have any other effect on the 
nutation. The theory and discussion by Woolard (1963) of the forced nutation is most 
complete. An excellent reference on this subject is Federov (1963). 

In 1876, Kelvin discussed the effects of earth deformation on the ocean and 
body tides. George Darwin and others attempted to measure the body tides. Darwin, 
in 1883, conceived of using the measured ocean tide itself, i.e., the difference between 
the absolute ocean tide and the absolute body tide, to establish the amount of body tide. 
Using long-period ocean tides (mostly monthly and semimonthly lunar tides), he found 
the amplitude to be two-thirds of the rigid earth tide and computed the rigidity of the 
earth of be the same as that of steel. A fine reference for body tides is Melchior (1966). 

The interaction of body tides, ocean tides, and polar motion has never been com- 
pletely understood; in some ways, our knowledge of the elastic and dissipative processes 
in the earth is not much advanced since the time of Kelvin and Darwin. Dissipation 
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presents a particularly formidable problem. The polar motion would decay to zero 
without continuous excitation, yet no energy source seems capable of providing sufficient 
energy (see Section 9). Similarly, body-tide observations with tidal gravimeters, 
though extensive, are not sufficient to indicate how much tidal lag there is in the solid 
earth. The time scale of the earth— moon system, given the current rate of lunar 
deceleration, is too small by a factor of 2 or 3. Dissipation from the solid earth 
would be continuous over geologic time, where the contribution from shallow seas 
would depend on their distribution. 

A general theory for ocean tides has never been formulated. The ocean tide is 
enormously complicated owing to irregular boundaries, loading, dissipation, and the 
necessary inclusion of inertial effects. Several recent papers are quite inconsistent 
in terms of the assumptions made and the results given. Compare, for example, 
Hendershott (1972) and Pekeris and Accad (1969); Hendershott and Munk (1970) give a 
general criticism of the subject. On the other hand, satellite orbits are sensitive to 
the ocean tide (see Section 8), and these tidal terms must be included to compute 
orbits of the highest precision. It is therefore imperative that the questions about 
ocean tides be resolved, perhaps by use of orbit analysis. 

We will see that satellite measurements allow some unification of the effects of 
the earth's elasticity. The determination of polar motion will be accomplished from 
kinematic effects, and the solid-earth and ocean tides will be determined from gravi- 
tational effects. 
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3. DYNAMICAL PROPERTIES OF POLAR MOTION 


The mathematical formulation of polar motion and body tides is most easily 
developed from the Liouville equation (Munk and MacDonald, 1960; hereafter referred 
to as MM, p. 10). Following MM, we adopt a set of geographic axes because they 
have a precise realization in terms of the position of observing stations. As long as 
the deformations and motions are small, we can treat the Liouville equation in the 
form: 


(3.1) 

( 

where 


i.+CT i = 0 d> 0 , 

1 r 2 r T 2 ’ 

V a r J l = - ff r*l > 

^3 = ^3 ’ 


„ _ C - A * _ d_ 

a -~ A ) dt ’ 


r A e 


(3.2) 


and 


♦l = “ 2 (“e c 13 + "e °23 + “e h l + h 2 ‘ L 2> ’ 

(C - A) O) \ 


<\>o - 


c 9<? - + w .h - h. + L. , 


v 2 .,2 \ e 23 e "13 w e 3 “1 "1 

(L - A ) w e 


(3.3) 


^3 = 2 ~ f "“e °33 “ “e h 3 + “e ) L 3 dt * ' 

“e C V 0 


The definitions are given in Chapter 6 of MM. 
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We have 


A 

l = 


'2 

1 

'3 


, the direction of the spin axis, 


or approximately the direction of the angular momentum of the earth; A, B, and C are 
the principal moments of inertia; c„ are the small contributions to the inertia tensor 
due to deformation; a^is the rotation rate of the earth; and h. and L. are the components 
of momentum and torque. 

In this study, we are not interested in the momentum and torque contributions, or 
the variation in the rotation rate. By introducing the complex variables 


i=i 1 + ii 2 , 


$ = 4^ + i(f>2 j 

where 

i = vrr , 

the equations reduce to 


(3.4) 


ii + or t = a cb , 
r r r 7 

a r^ = A [ w e C 13 + C 23 + l(w e C 23 " C 13 ) ] 

• * 

Assuming a rigid body, we have c ^ = c 2 3 = c 2 3 = c i 3 = 0 


(3.5) 


WV1 

i - iff i = 0 , 

r ’ 
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which has the solution 


i = c x e + c 2 ; 

c, c are constants. The free period 
J. ^ 

r= A/(C - A)/co e .... . 

If co g ~ 1 revolution per day and (C - A)/A is obtained from the precessional constant 
H = (C - A)/C = 0. 0032729, i. e. , (C - A)/A = (1 - H)(C - A)/C = 0. 003262, we have 
a free period of 306. 5 days, the result Euler obtained. 


The deformation of the earth is a complicated problem in elasticity. We consider 
forces as derived from disturbing potentials. A. E. H. Love introduced some dimen- 
sionless parameters that summarize some of the earth's elastic properties. Their 

detailed theoretical computation has been carried out by Takeuchi (1950), Longman 

* 

(1963), and Farrell (1972). These numbers are defined in the following way. The 
displacement of the ground due to a potential of degree n is 


n (surface) 
g n 


and the additional gravitational potential due to the motion of mass is 


k U 


n n 



2n+l 


(3.6) 


(3.7) 


If the force loads the earth, e. g. , the ocean, then h^ and k^ are used. The response 
of the earth to deformation certainly changes with position. Kaula (1969) considers 
the Love numbers expressed in spherical harmonics. 


There is no direct way to measure k and h by terrestrial means. Only the linear 
combinations 1 + k - h from the deflection of the vertical and 1 + h - (3/2)k from 
gravity. However, k£ is determined from the free nutation. The results are not 
altogether consistent, and interpretation is complicated by the frequency dependence 
of the Love numbers owing to resonance effects with the core. We can say that 
k 2 «0.3 and h 2 =0.61 and that k < 0. 504 h 2 (Melchior, 1966). 
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The displacements due to rotation about 9 , 9 , 9 are the same as those due to 

1 Z O 

the potential of second degree: 


1 2 . 2 , 2 , 2 , , , 

2 w g ( x + y ) - co e z(i lX + ^ 2 y) , 


x, y, and z being the coordinates of the displaced point. The first term adds to the 
oblateness and does not concern us. The second term gives rise to the additional 
external gravitational potential 


-k 2 go 2 z(i jX + i 2 y) 



(3.8) 


MacCullagh's formula (Jeffreys and Jeffreys, 1956) gives the gravitational 
potential due to a deformed earth: 


V = G 


(A + B + C) r 2 - 3(Ax 2 + By 2 + Cz 2 + 2c^yz + 2 c„xz + 2 c„ X y 


'13 


'12 


2 r‘ 


(3.9) 


Comparing terms, we have 

2 a 5 

w a 

p = k e e o - a 9 
C 23 k 2 3G 2 ai 2 ’ 

gj 2 a 5 

C 13 = k 2 3G 1 1 = ai l ’ 


(3. 10) 


a =.k, 


2 o 5 

gj a 
e e 

2 3G 


Inserting these into Eq. (3. 5), we have 


O’ 


a r <|} “A [w e i l + i 2 + i(w e i 2~ i l )1 


(7 

= A ^ u e J " 


(3. 11) 
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giving 


i(A + a) i + u (C - A - a) l = 0 . 

This differential equation has the solution 


(3. 12) 


r. c, e-“/ T + e 2 

with the free period 

A + a 

T ~ OJ (C - A - a) 
e 


(3. 13) 


(3. 14) 


If we now take the free period as given by latitude observations, then the Love 
number is determined to be k = 0. 29. Further consideration of these equations could 
shed more light on the dynamics of polar motion. 


Before going on to the main body of the report, let us raise several points that 
will be useful later. Kelvin has shown (MM, p. 29) that for an incompressible 
homogeneous sphere of rigidity p, 


k = | —~2 , (3- 15) 

0 (n - 1){1 + p[2(2n + 4n + 3)/l9n]} 

where the dimensionless rigidity p is related to the rigidity by 


V- 


19 

2 


H— * 
Pg & e 


2.3 


(3. 16) 


It is from this formula, by use of the Love number obtained from the Chandler motion, 
that the rigidity of the earth was initially determined. 


We can introduce dissipation into these equations by use of a complex Love num- 
ber (MM, p. 153): 
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U singly in the previous equations includes the effect of dissipation. 

It is now instructive to consider a simple solution of these modified 
equations. If we let 

2 5 




co a 
e e 

CT = 

Y 

3G 



2 5 
co a 

K = 


e e 

k 2 

3G 

r = 

i±_ 



2Q 

7 

A' = 

= A + K , 





where 


) denotes the real part of () , 

then Eq. (3. 12) takes the form 

(A' - iTK)T- i(C - A' + iTK)7= f(t) . 


The homogeneous solution is, of course, 


i O t 
c 


H - (const) e + another constant , 


where 


C - A' + iTK 
c A' - iTK 


a = i} a = A'(C - A') - 1K D 
n JiU c (A' - iKT)(A' + iKJ 


Kr) 


(3. 17) 


Liouville 


(3. 18) 


(3. 19) 


(3.20) 
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Now, for a forcing function 


ia,t 

f(t) = a e a 


we have the particular solution 

ia, t 

T = ia_e 

(A'- iKT)[a d - or +i(Kr/A')] ' 

We can write the real part of the frequency as 


(3.21) 


a = [(C - A')/A'\ - (FK/A') 2 _ (3 ; 

n [1 + (TK/A') 2 ] 

For k 0 = 0. 3, cj = 0. 7X 10" 5 sec" 1 , a = 6. 3 X 10 8 cm, G = 6. 67 X 10" 8 , p = 2. 3, 
4 e aa e _q 

Q = 30 to 500, A~B~C~8X 10 , and (C - A')/A' ~ 1/430 « 2. 3 X 10 , we have 

-5 

FK/A' « 10 . Therefore, the dissipation has no effect on the free period, and 

further, the amplitude of the forced motion is controlled by the KT/A' term. 


The natural decay time is also determined by 


rK „ k 2 ^e a e _ k 2 ^e a e ^ 

A 6 G Q (C/Ma 2 ) Ma^ 6GQ(C/Ma^) 

1.2 X 10~ 3 

Q 

The decay time is 


T = x i° 3 toys = 5-^5 years . 


(3.23) 


(3.24) 
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Finally, we wish to make a correspondence between the contribution of the defor- 
mation to the external gravity field expressed by Eq. (3. 9) and the form of the poten- 
tial used to develop orbital perturbations. It is customary to write the external 
potential as 


V = 


GM 


°° l JL 

5 <m (sln <•>> cos mX + sin mX > 

L 1=2 m=0 


(3.25) 


where the P^ m (sin cf>) are fully normalized associated Legendre polynomials, i. e. 


/ p L (sin « 

Sphere 



mVcos dcj) d\ = 4tt . 


(3. 26) 


If we let x = r cos <|> cos \, y = r cos <f> sin and z = sin <j> in Eq. (3. 9) and compare them 
with the expressions in terms of Legendre polynomials, we find, for example, 


P 21 (sin <t>) = 3 



sin 4> cos <(> , 


yz = 


2 


r 


2 _ 

cos 4> sin 4> sin \ = ■ P . (sin c}>) sin X . 

yjl5 


Therefore, we have 


C 21 = 


S 21 


i 2 3. 

- k„ oo a i. 

2 e e 1 

y[T5 GM 
, 2 3. 

“VeV|2 

•\/Ti GM 


(3.27) 
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4. POLAR MOTION AND EARTH SATELLITES (KINEMATICS) 


We now turn to the observable effects of polar motion on satellite orbits and 
satellite observations. We will consider first the kinematic effect on observations 
of satellite positions. Then we will explore the orbital perturbations arising from 
variations in the mass distribution, i.e., gravitational perturbations. 

Consider a reference system established with respect to some arbitrary point 
on the earth's surface. This point is chosen to be near the spin axis of the earth, 
and we will choose for this discussion the Conventional International Origin (CIO) . 

Aside from tides and crustal motions, which we ignore here, the observing stations 
are rigidly fixed with respect to this point. Each station has its geographic latitude cj>. 
The realization of the position of several stations in this reference system is the 
subject of geodetic research. Satellite -tracking data have been very important in 
achieving a consistent set of coordinates (Gaposchkin and Lambeck, 1971). The actual 
values of the station coordinates establish the origin in an analogous way to the mean 
observatory of the International Latitude Service. 

Consider one such station, point A in Figure 1. The terrestrial system has its 

A 

pole (the CIO, for example) in the direction of P, which is at an angle x with the spin 

axis <£). The spin axis is approximately the same direction as the angular momentum 

L of the earth. It is the angular momentum of the earth that is constant, as we are 

assuming no torques. In the classical Poinsot construction, the spin axis is f • x 

from the direction of L, where f is the dynamic flattening. With l/f = 298, and 
-6 -9 

x~0.5X 10 radians, the difference is 1. 6 X 10 radians or about 1 cm at the 
pole. We assume the spin axis is fixed in space. 

Consider a satellite with inclination I. The inclination is defined with respect to 
the inertial reference system defined by L. If we assume that the satellite's orbit is 
well known and that it can be observed at any point in the orbit, then we make an obser- 
vation from A to S and 12 hours later from A' to S'. In the first case, the angle AOS is 
cf> + I - x, and in the second, <f> - I - x. Using the adopted values of <j>, we determine x. 
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This is the projection of the pole position on the meridian of A. A station 90° east or 
west would determine the other component. We note that cj> - x is an apparent change 
in latitude that leads to the occasional use of latitude variation as a synonym for 
polar motion. 

A polar satellite would be used in the same way. Instead, we would determine 
v + oj - x + I and v + to - x - I. We see that a satellite at any inclination can be used, 
the accuracy depending on that of the combination of along-track v + w and cross -track 
I components. The inclination of a satellite is inherently more accurately known, so 
the optimum procedure is to make observations at the point highest in the orbit when 
this point is in the meridian of the station. 

This analysis has been simplified by considering only observations in the meridian. 
Generally, a station will have observations through 10 to 20° longitude, and a single 
station will provide some information about both components. This will be more 
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apparent when we consider determination of pole position as part of the process of 
orbit determination. 


Finally, since the polar motion gives rise to diurnal variations in the apparent 
latitude <j> - x, the principal corruption in pole position determination comes from 
diurnal errors. There are important diurnal perturbations in the satellite position 
from the first-degree tesseral harmonics. Therefore, these H = 3, 4, 

geopotential coefficients must be especially well determined. 

Ultimately, when orbital uncertainties are reduced, the possibility of diurnal 
effects in the observations will become important. Environmental conditions at a 
station have very strong diurnal variations. 


To determine the pole position in the context of orbit determination, we consider 
the expression relating the observation and the computed satellite position 


p = r - R 


(4.1) 


Here, r is the position of the satellite and R is the position of the observing station. 
The observation equation for direction or range observations is derived from a projec 
tion of this expression. For example, the range is determined from 


A A < n\ 

P = P • P = P • (r - R) 


(4.2) 


The determination of any parameter p entering r or R is obtained from 9p/9p. The 
coordinates of the pole x, y enter R as 



cos 0 -sin 0 0 


r~ 

X 

i 

o 

«— l 

R = 

sin 0 cos 0 0 


0 1 y 


_ .0 0 1_ 


_x -y l_ 


X , 


(4.3) 


where 0 is the sidereal angle, the X are the rectangular coordinates of the observing 
station in an earth-fixed system, and the x, y are the coordinates of the pole given in 
the convention of the Bureau International de l'Heure (BIH). The coordinate x is 
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positive along the Greenwich meridian and positive y is to the west, i. e. , in the 
opposite sense to the y coordinate of the station. The partial derivatives easily 
follow: 


i£ _ 
ax 


“z cos 0“ 


z sin 0 

z sin 0 

i£ _ 

’ 3y 

-z cos 0 

-x _ 

_ y 


(4.4) 


It remains to consider if these observation equations are in themselves singular 
and, further, if the resulting set is singular when they are augmented with equations 
for orbital elements. Both questions can be answered negatively from the following 
cons ide rations . 

The normal equation for the system itself for one station is 


n 


n 


y, (z 2 cos 2 0. + z 2 sin 2 0. + x 2 ) ^ (z 


2 2 
cos 0. sin 0. - z cos 0. sin 0. 
ii li 


n 


n 


Z 2 2 V > 2 2 2 2 2 

(z cos 0. sin 0. - z cos 0. sin 0. - xy) y (z cos 0. + z sin 0. + y ) 


n 


2 2 
z + x 


L -xy 


-xy 

2 , 2 , 

z + y J 


(4.5) 


with determinate 


TJ 2 2 

nR z , 


which is zero only if the station is on the equator. This does not mean that we cannot 
determine polar motion with an equatorial station, only that we cannot separate the 
x and y components. The component along the meridian can be determined. 
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We observe that the correlation between x and y is proportional to xy and can be 
reduced to zero by placing the station with either x or y = 0. This is intuitively 
reasonable, as the component along the meridian of the station is most easily deter- 
mined. By choosing x or y = 0, one decreases the correlation and improves the 
determination of the component along the meridian at the expense of the accuracy of 
the orthogonal component. The element of the normal system corresponding to the 
orthogonal component is minimum and results in a maximum variance in the inverse. 

Since the data are taken at arbitrary times, the geometry of the observations 
controls the part of the normal system relating the orbital elements. In general, we 
can expect (and in fact find) no singularity in determining pole position with orbital 
elements . 
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5. POLAR MOTION AND EARTH SATELLITES (DYNAMICS) 


In the theory of free nutation, we found that the elastic deformation that lengthens 
the period of free nutation also gives rise to tesseral harmonics C ^ and S 21 , which 
depend on the pole positions l ^ and l Can we hope to determine these coefficients 
by satellite perturbation analysis and recover the pole position? 


The orbital system defined by the direction of L is the system in which we study 
orbital perturbations. The gravity field is expressed in the system of principal axis. 
We want to explore the consequences of orbit determination in which the geopotential 
is moving. The dominant part of the anomalous gravity field is due to J 2 = -/5X 
We assume the system of principal axis has coordinates £, 77 with respect to L. In 
the principal-axis system centered on |, 77 , the gravity field can be approximated 

by 



1 + c 


20 



P 20 (Sin ^ 


(5.1) 


We could use the general transformation of spherical harmonics to find the expression 
for V in the system oriented along L. However, if we assume £, 77 to be small, a 
simple expansion suffices and we have 

P 20 (sin 4»') = /5 sin 2 4 ,' - | j = V5 (§ z' 2 " • 

In the orbital system, 

z' = z - |x - Tjy 
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and 


3 2 1 


P 2o( sin ^ = \ 9 z - 9 - 3^xz - 3rjyz + 


2 2 


= P^sin 4>) - £>/3 P 01 (sin cj>) cos X - rjV’3P gi (sm $)■ sin X , 


21 


21 ' 


(5.2) 


so we have 


V = 


GM 


1 + C 


20 


0 )' [ 


P 2 Q(sin <j>) - /3 P 2 ^(sin <)>)(£ cos X + 17 sin 


in X)] 


(5.3) 


Therefore, the motion of the principal axis appears as tesseral harmonics of 
degree 2 and order 1: 


C 21 C 20 ^~ 3 ^ ’ 

(5.4) 

®21 = ' °20 V ’ 3 * • 


The gravitational potential sensed by the satellite in orbit referred to L contains the 
terms 




/ 3 | - 


k u 2 a 3 f 
2 e e 

/I5 GM 




C 20 /3 * 


i 2 3, 
k„u a St 
2 e e 

'/l5 GM 


2 


(5.5) 


Using quantities known with sufficient accuracy, we have 


C 21 - 0.838X 10~ 3 £ - k 2 0.893X 10 3 i , 

S 21 = 0.838 X 10“ 3 rj - k 2 0.893 X 10~ 3 i 2 . 


(5.6) 
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If we can choose the origin of the geodetic system to be the same as the principal axis, 
which is not the case for the CIO, then £ = i rj = f 2 , and the values become 


C £1 = (0.838 - k 2 0.893) l X 10~ 3 , 


S 21 = (0. 838 - k 2 0. 893) i X 10 -3 , 


(5.7) 


with the elasticity reducing the effect by about one-third. The minimum effect would be 
zero for k = 0. 936, i. e. , if the Love number was the fluid Love number (MM, p. 26). 

-5 

If we assume i^=| = 0.5X10 and k 2 = 0. 30, we have 
Cgi = -3 X 10~ 9 . 


Until now, this accuracy has been achieved only for the zonal harmonics of the geo- 
potential (Kozai, 1969). Zonal harmonics cause long-period and secular perturbations 
that can be averaged for months or years. For close-earth satellites, the tesseral 
harmonics and give rise to diurnal perturbations. These perturbations will 
have the same frequency spectrum as do all the first-order tesseral harmonics. 
Synchronous satellites are also sensitive to the first-order coefficients. These 
synchronous satellites will be in the resonant state and the orbital effects considerably 
amplified. The use of synchronous satellites to study the polar motion has been pre- 
viously discussed by Gaposchkin (1968); 


We consider the use of laser retroreflector satellites as the optimum method for 
determining the geopotential today. Table I contains the amplitude of the perturbation 
for the seven existing laser reflector satellites and the nominal orbit of Geos C. The 
perturbation is estimated by using the approach given by Gaposchkin (1970). 


The principal perturbation due to and has a diurnal variation, which is 

the same frequency as the kinematic effects (see Section 4). The 1-m perturbation 

“5 6 

comes from a 0.5X 10 0 X 6 X 10° m = 30-m polar motion. We can imagine a process 
where we determine i ^ and i from the kinematic effects and then use them in 
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Eqs. (5. 7) with a nominal value for k . We can then redetermine SL and i , using 

„ ^ J. ^ 

these values of Cg^ and Sg^. Since the error in computing the perturbation is l/30 
the error in the assumed position of the pole, the process will converge. 


Table 1. Perturbations due to C ^ and Sg^ with l = 0. 5 X 10 5 and k = 0. 30. 


Satellite 

Inclination 

a (Mm) 

Eccentricity 

Amplitude (m) 

6406401 

BE-B 

79° 69 

7. 360 

0.013 

0.70 

6503201 

BE-C 

41. 17 

7.502 

0.025 

1.41 

6508901 

Geos 1 

59.36 

8.074 

0.071 

1.07 

6701101 

D1C 

39.96 

7.319 

0.051 

1.51 

6701401 

DID 

39.43 

7.603 

0.084 

1.40 

6800201 

Geos 2 

105. 80 

7.708 

0.032 

0.88 

7010901 

Peole 

15.0 

6. 999 

0.017 

1.20 


Geos C 

67 

7, 28 

0.006 

1. 12 


The analysis of pole position, by using values of and determined by satellite 
perturbation analysis, can be approached as follows. We imagine employing a satellite 
with minimum sensitivity to the gravitational perturbations, i.e., BE-B or a high, 
heavy satellite such as Lageos, to determine the pole position. Then, by using the 
gravitational perturbations of a sensitive satellite, such as BE-C, D1C, or DID, we 
can determine Cg^ and Sg^ as functions of time. We can then: 1) Assume a value of 
k 2 and determine £, r) from Eq. (5. 5) or (5.6). 2) Assume a value of |, rj, say 
i 1 = |, l 2 = n> and determine kg. The temporal variation of kg could be investigated, 
(Gaposchkin, 1968, 1972)* 3) Assume a particular form for the variation of £, r/, e.g., 

4 = cos oot, n = t]q sin cot, and determine £q, t/q and kg. This would require a 
long series of data, several oscillations, i. e. , several years. 
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6. BODY TIDES AND EARTH SATELLITES (KINEMATICS) 

The application of an external potential of degree n, U n (r,cf>, X), causes the surface 
of the earth to move vertically by an amount 


Ar 

n 


h 

= fV n (r,4>,\) 


(see Section 3). Now, U 2 can be written (MM, p. 68) as 

a 4 

U 2 = g 2 "M ^3 b £ «>> 008 13 ■ 

The largest term is M 2 with b = 0. 908, and we have 

2 

Ar M = h 2 (53. 7 cm) 0. 908 £2^L_3L C os 2 (X Q - - 79° 8 + X) . 

2 

By using the rough value of h = 0.61, 

z 

Ar^ = 14. 9 (cm) cos^ cj) cos 2 (\q - X' - 79° 8 + X) 

2 

We see that the kinematic semidiurnal raising and lowering of the observing station 
is 14. 9 cm, which is maximum on the equator. Takeuchi (1950) has computed the 
complete deformation by use of a model earth and obtains 22 cm for the total change. 
We conclude that only when observations and data-analysis capability reach centimeter 
accuracy will this effect be observable. The current value for h 2 has sufficient accu- 
racy to be used as a preliminary number in the reduction of observations. 


23 



The previous calculation was based on uniform properties of the earth. There 
are local variations; for example, tectonic conditions can control the deformation 
due to atmospheric pressure. With certain pressure distributions, the vertical 
movement amounts to 15 cm in northern Europe (Tomaschek). 
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7. BODY TIDES AND EARTH SATELLITES (DYNAMICS) 


We have defined the additional potential due to the deformation from a potential 


U as 
n 


AV = k U — 
n n n \ r 


2n+l 


(7.1) 


The total potential acting on the satellite is then 


1 + k — 
n \ r 


2n+l 


U 


n 


(7.2) 


The potential that acts on the satellite owing to the moon (or sun) can be written 



(7.3) 


where r and r' are the position of the disturbed body (the satellite) and the disturbing 
body (the moon or sun), respectively. We can write 


1 _ 

A 


oo a 



i=0 m=0 


1 

21 + 1 


ITT P (m (s “ « P im (sln *'> e im<W> 


(7.4) 


where r, cf>, X are the spherical coordinates of the satellite and the primed symbols are 
those of the moon (or sun). To calculate orbital perturbations, we will use the 
gradient of Eq. (7. 4) with respect to the satellite position, and we can drop the i = 0 
term. The i = 1 term just cancels the (r • r')/| r' | , so we have for the third-body 
potential 
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(7.5) 


V 


ou ~ 

= 2iVr 


1=2 m=0 


i 

■7+T P im (sin +> P im (s “ +'> 
r' 


e im(\-\') 


Therefore, the potential acting on the satellite, including the tide, is 


00 n 

t'TPW 1 


k, a 2,+1 
i e 


V'- GM 'Jlr L, 2s 2 TTT l 3ITI + 3T + i>i ; p im< sin « p im< sin +'> e 
i=2 m=0 


n im(\-\') 


(7.6) 


To include the effects of the tidal phase lag, one introduces a fictitious moon lagging 
the real moon by At. In this case, we cannot write the disturbing potential in such a 
compact form. We proceed, assuming At = 0, i. e. , no phase lag, the revision of the 
theory being relatively straightforward if lag is desired, either by use of At or a 
complex k n> 

To develop perturbations in satellite position, we have to express the potential 
in orbital elements. The first step is through the relation (see Appendix I) 


= . . ,, imX 

P Jm (sm e 




p=o 


D (i) e^-ZpHv+wJ+mn] 
imp' 


•which gives 


(7.7) 


V'= GM' 


00 i t t 

a/ +m 1 

f r‘ , k ^?:M 


2i + 1 1 

l r ,i+l r ,l+l r l+l/ 


°imp® D i(-m)p 


,(!') , 


(7.8) 
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where 


ip= (£ - 2p)(v + to) + (i + 2p')(v' + a/) + m(ft - Q') , 


by noting 

Vm)^ * 

The final step in expressing Eq. (7. 6) in Kepler elements involves expressing 
rPQmv as a series in a; Kaula (1966) gives formulas for 


1 i[(i-2p)v+$ 

T+I e 

r 


- 1 V r m P i[(f-2p+q)m+^] 

“a^+l 2~i W 


(7.9) 


A more general form is given by Plummer (1918, p. 44): 


n imv 
r e 




iqM 


q=-°o 


(7. 10) 


where x"* n 1 (e) are known as Hansen coefficients. 

c[ 

for X^’ (e) given in Plummer, where 


There is a classical development 


x"’ m =(l + |3 2 )Vj <qe)x"< m , 
q Z-i p q,p ’ 


j, m . a . q-p-m /n + 1 - m\ _ . , . , , _2, 

q’ p = (_(3) (q-p- m J F(q ~P" n ‘ 1 > _m ~ n ~ 1 > q " P " m l » p ) » 

q - p - m > 0 , (7. 11) 

(Eq. cont. on next page) 
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vii, m . Q . -q+p+m /n + 1 + m \ „ , . . , , , . „2. 

qip =( " W U + p+mj F(-q + p-n-l, m - n - 1, -q + p + m + X, |3) , 


q - p - m < 0 , 


- n -1, -m - n -1, 1, (3 2 ) , 


q - p - m = 0 , 


(7. 11) 




in the notation for a hypergeometric series. These may suffice for theoretical 
studies. For practical problems, we prefer the expression as a polynomial in e. 
Andoyer (1903) developed a recursive method for computing this polynomial, which was 
rediscovered by Izsak et al. (1964). Chemiack (1972) gives the practical details 
for computing the X^’ m (e) . 

By using Eq. (7. 10), V' becomes 


oo £ 9 £ oo oo 

-ZELL £ £ V £mpp'qq' ’ 

9.-2 m=0 p=0 p'=0 q=-oo q'=-oo 

where 

(-l/ +m 

V £mpp'qq' " 29 + 1 D £mp^ D £(-m)p^ 


X 


2f+l 


.£+1 q 


m (e) X'f _1 ’ m (e') 4* —~~t 


k 9 a e v -£-l,m . v -£-l,m. , 
,£4-1 J 4-1 q (e) q' (e > 


e l * , 


a' a 


and 


i/j = qM 4- q' M' 4- (£ - 2p) oj 4- (£ - 2p') J 4- m (S2 - n') 


(7. 12) 



Considering the size of the perturbations, we chose to use a linear theory 
developed in terms of the Lagrange Planetary Equations (LPE). These are given in 
Kaula (1966): 


dto - cos I 9R 

dt “ 2.. 2, 1/2 . T 91 

na (1 - e ) ' sm I 
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H - e ¥/ 2 

2 

na e 


9R 

9e 


> 


d£ _ 1 9R 

dt 2~. 2,1/2 . j 91 

na (1 - e ) ' sm I 


dl _ cos I 9R 

dt 2~ 2, 1/2 . T 9 co 

na (1 - e ) ' sm I 


1 

2_i 2,1/2 . T 

na (1 - e ) sin I 


9R 

9Q 5 


de 

dt 


2 

na e 


9R_ (1 - e 2 ) 1//2 

9M 2 

na e 


9R 

9co 


> 


dM _ 1 - e 2 9R 2 3R 

dt n 2 9e na 9a 
na e 


(7. 13) 


da _ _2_ 9R 
dt na 9M 


By using Eq. (7. 12) for R in (7. 13) and integrating, we obtain the lunar or solar 
perturbations. They have been developed by Gaposchkin et al . (1973) by use of computer 
algebra, and we will not give the full expression here. From the argument ip, we 
observe that there is a rich spectrum of perturbations. The orbital effects are large 
for those terms involving slow rates, i. e., q = 0. 

With respect to tidal studies, Kozai (1970) has used the £2 - £2' terms, which have 
a very long period. Using these terms requires careful study of all long-period 
effects. The principal difficulties arise from radiation pressure and drag and are 
complicated by our imperfect knowledge of the satellite aspect, characteristics, and 
the atmospheric density. For this reason, we concentrate on the two-week lunar 
tides. Though the orbital effects are smaller, the other effects are manageable and 
the theoretical difficulties are reduced. 
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We estimate the size of the direct lunar perturbation as 


« 1. 5 X 10 5 « 120 

2nn' 


m 


where 


M' 

^ M ’ 

which is quite a good estimate. The tidal effect is k 2 (a J 2) times the direct effect 
and is about 15 m. 
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8. OCEAN TIDES AND EARTH SATELLITES (DYNAMICS) 


The observations of ocean tides are extensive, but mainly along coastlines. The 
theoretical treatment of ocean tides remains a fundamental problem. Hendershott and 
Munk (1970) give an excellent review of this subject. For our purposes, we note that 
the bulk of the tide observations are taken on coastlines that have very important 
boundary effects, and that all treatments existing today are inadequate in one way or 
another. Although there is superficial agreement in the cotidal charts, the spherical 
harmonic coefficients differ by more than 100% and are often different in sign. Only 
a few terms in the general tidal development are important for satellite perturbation 
analysis, since in exactly the same way that the satellite is a frequency filter for 
the gravity anomalies, it selects only a few of the tidal terms. 

Most of the tidal theory has been worked out for the M tide, the principal semi- 
diurnal tide (Pekeris and Accad, 1969; Hendershott and Munk, 1970; Hendershott, 

1972). It is the two-week tide, in the inertial reference system, that would arise from 
the moon in the earth's equator. 

Current determinations of the low degree and low-order terms in the tide are quite 

inconsistent. An equilibrium theory, which surely is inadequate, predicts the 
i2\ 

P (sin4>) e tide of 30 cm. Pekeris and Accad give a solution in terms of a number 
map. Lambeck and Cazenave (1973) have determined by harmonic analysis an 
amplitude of 4.4 cm and a phase of 330°. Hendershott gives an amplitude of 50 cm. 

We can only comment that these estimates roughly agree as to the order of magnitude. 

For further analysis, we assume a tide raised by a disturbing body in a reference 
system centered at the earth's center and oriented toward the disturbing body. We 
have l/Afrom Eq. (7.4), and the tide raised is developed in spherical harmonics: 





(sin <j>") 


imk 

e 


n 


9 


( 8 . 1 ) 
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where 


?im C im" lS im ‘ (8,: 

If the tide is static in this system, it follows the disturbing body and fim is 
constant. If there is no lag, then = 0. The external potential due to this tide, 
including the loading effect, is 



(l + k i)4^p u a' +2 — 


(2i + 1) r 


i+1 




lm P im (sl1 ♦"> e 


imV 


(8.3) 


where p is the density of the water, k', is the Love number for a deformation that 

CO 

loads the earth, and and X" are the coordinates of a point referred to the position 
of the disturbing body. The quantity <j>" is measured from the lunar orbit and X" is 
measured from the position of the moon v' + u/ in the lunar orbit. If X' is referred to 
the lunar node, then X' = X" - v' - co' . We can express this potential in the inertial 
system by using the transformation (see Appendix 1) 


Stm 


(sin <(/) e 


i[mX'-m(v'W)] 


<i) S - m E fms H') P |s (sin 4.) e «^- m <G'^'>] 


s =-i 


( 8 . 4 ) 


We can now express the potential due to ocean tides in the orbital elements of the 
satellite as 


1=2 m=0 s=-i p=0 

where 

v _ r m i+s /n m F i-n n m e i[(i_2p)(v+w)_m ( v ' +w ' )+sfi ~ mSri 

V imsp ( 1} ^Ims' L > U fsp W e •» 

( 8 . 5 ) 
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in which 


r im 


^Cpjl + k') a 


t+2 


2i + 1 




4+1 r im 


( 8 . 6 ) 


Similarly, if we can take the definition of the M 2 tide with I' = 0, we have 



4= 2 m=0 p=0 


where 


m . i[(-^-2p) (v+^+m^-v'-co'-n')] 

V lmp =r (m li) ' “imp ® 6 


(8.7) 


With this expression for the potential, we can use the Hansen coefficients, Eq. (7. 11), 
to express it in Kepler elements and then put it into the LPE, Eqs. (7. 13). Inte- 
grating the resulting equations gives the perturbations due to the ocean tides: 



i + 1 4 ttG p (1 + k') 

± OJ j- 


4+2 21+1 
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e 



( 8 . 8 ) 


ip = qM + q' M' + (4 - 2p) co + m (£2 - £2' - a/) , 

ip = qn + q' n' + (4 - 2p) oo + m (h - fr - w' ) , 


(8.9) 
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( 8 . 10 ) 

(Eq. cont. on next page) 
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We have given the formulas for the M tide caused by a fictitious moon in an 

Li 

equatorial orbit. These formulas are easily extended to a tide oriented in the plane 
of the lunar. orbit. We would use the potential given by Eq. (8. 5) . We acquire another 
summation index, and the formal expressions for the perturbations would be 

6£ impqq's = 6£ £mpqq' E ims^ _I ^ ’ * 8 ' ^ 


with 


tp = qM + q' M' + (S. - 2p) u>+ sf2 - m(co' + Q') 


and 


ip = qn + q'n' + (i - 2p) <L + sf2 - m (w' + f2') , 


with £ a generic element. 


It has been convenient to develop the theory for orbital perturbations in complex 
notation. It is easy to utilize these formulas as given by automatic use of complex 
variables in modem Fortran compilers. To work in terms of real variables, one 
can easily separate the real and the complex parts of ^ m , (i) , k^, e the other 

variables that appear are real. 


It is instructive to determine the ocean-tide equivalent of the body tide; however, 
we can do this only approximately. We make the correspondence by comparing the 
potentials in Eqs. (7. 8) and (8. 7) for a particular imp combination. We have 
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and 
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where 


(8. 13) 


ip = (£ - 2p)(v + oo) - m(v' + cj' + n' - G) . 


We note that the lunar inclination is I' - 23° ± 5°, and that = 0- 925, 

D 2(-2)l = P 2(-2)2 = 0 *®036. So f° r principal semidiurnal term, we can 

take £ = 2, m = 2, £ - 2p = 2, p = 0, and p' = 0, which gives 


4irGp 
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or 
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co 


where k 2 would have a complex value as defined by Eq. (3. 17). Using typical values, 
we have 


k 2 = 0 - 0114 £n/ D 2(-2>0< I '> 


(8. 16) 


From Lambeck and Cazenave (1973), the Pekeris and Accad (1969) solution with 
dissipation gives 



-i 330 it/ 180 

4*4 © 


= -2. 19 - i 3. 81 (cm) . 


We then have k = -0. 026 - 0. 047 i. If we then add this to the body tide, we 
obtain the effective Love number that a satellite would sense. Choosing 
jbody _ q_ go no dissipation, we have 


36 



.effective _ . body .ocean 
z _ k 2 + k 2 

= 0. 274 - 0.047 i . 

We would conclude that the Love number was 0. 276 with a phase lag of 9° 73, or 
38. 9 min. 

Conversely, we could determine the ocean tide if we were to adopt a value of the 
effective Love number and the body Love number. 

We have analyzed the P component with a 14-day period. The ocean tide has 
a much richer spectrum than does the body tide, with important contributions from 

and Pg 2 - They will give rise to similar perturbations to the P^, in both ampli- 
tude and frequency. Therefore, the oceanic contribution to the satellite perturbations 
must include these terms, and we cannot easily determine just one 

Analysis of tides by Anderle (1971) and Smith, Kolenkiewicz, and Dunn (1972) used 
these two-week terms. Newton (1968) and Kozai (1970) used the terms in £7 - £7'. The 
perturbations due to those terms from the ocean tide are much smaller, owing to the 
expansion in Hansen coefficients, and should not influence the determination. There- 
fore, we can roughly reconcile the determinations of Anderle, Kozai, Smith, 
Kolenkiewicz, and Dunn. Newton's result could be explained by the fact that he used 
^2 = ^2’ w ^ ereas ^2 “ ~^2’ em phasizes the importance of atmospheric tides in 

the treatment of long-period perturbations. 
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9. ANALYSIS OF EXISTING POLAR-MOTION DATA 


The availability of information on pole position from 1846 through 1970 invited 
a rediscussion of the data and properties -of the polar motion. The results of the 
analysis have been reported by Gaposchkin (1968, 1972). 

Since completing the analysis, I obtained from Madame Rykhlova the computed 
pole positions published in graphic form, which we read to obtain tabulated values. 
Reprocessing these improved data gave better results in the sense that the high- 
frequency-noise part of the spectrum was very much reduced in amplitude. The 
properties of the motion and the amplitudes remained within the expected uncertain- 
ties. 


Since the analysis was completed, the complicated spectrum in the Chandler Band 
has been studied in other ways. The polar motion causes an oceanic tide called the 
pole tide. This pole tide has been studied with a frequency spectrum analysis, and 
the two main frequencies have been found. Clearly, the pole tide, being excited by 
the polar motion, would exhibit the same spectrum. However, this corollary 
evidence rules out the possibility that the complicated frequency spectrum is in some 
way due to some property of the reduction methods. One such study has been carried 
out by Smets (1970). 



10. ANALYSIS OF SATELLITE -TRACKING DATA 


The 1969 Smithsonian Standard Earth (II) was determined with a combination of 
camera observations, simultaneous tri angulation, and dynamical data, from surface- 
gravity data compiled in 1964. hi addition, the Deep Space Network of the Jet Propul- 
sion Laboratory provided observation equations based on some deep-space-probe data. 
A limited amount of laser data were also used. 

0 

The camera data had an assumed accuracy of 3 to 4 arcsec. At 2 X 10 m, this 
corresponds to a 30-m observation. The station coordinates were determined to 
between 5 and 10 m. This figure was given at the time and has been verified subse- 
quently by comparison with the Smithsonian Standard Earth (III) . The orbit computation 
seemed to be ~10 m for the Geos 1 type satellites, though this proved to be optimistic. 
When sufficient laser data for good orbital coverage became available from the ISAGEX 
program, the orbital fits were 20 to 30 m. One significant contribution to this uncer- 
tainty is the tide perturbations. Though the tides had been studied by Kozai with long- 
period terms, no short -period, 14-day, tidal terms were used in the analytical theory. 
Only with the complete treatment of lunar and solar perturbations by use of computer 
algebra have we satisfactorily accounted for the body tides. In this report, we give 
the necessary formulas to include the effects due to the body and ocean tides. 

The attempt to determine polar motion in 1969—70 by use primarily of photographic 

data was unsuccessful. The investigations were, successful, however, in exposing 

*..■» •. .•/' '-I i %f ■ s 

several aspects of the analytical" process that required refinement. With the advances 
that have been made in orbit computations through a revised gravity field, improved 
station coordinates, and a complete tidal treatment, and with the improvements that 
have been obtained in the accuracy, reliability, and coverage in significant bodies 
of laser-ranging data, we can expect that future efforts will produce useful deter- 
minations of polar motion. 
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APPENDIX 1 


ROTATION OF SPHERICAL HARMONICS 

Legendre functions arise naturally in the solution of Laplace's equation in 
spherical coordinates. They are also used as a complete set of orthogonal base 
functions for mapping arbitrary functions in spherical coordinates, particularly in 
potential theory. As a set of base functions, we naturally need their expression, 
transformed from one coordinate system to another. We consider here the rotation 
of a coordinate system. This transformation has been used in quantum mechanics, 
and proofs of the result described here can be found in that literature. For our 
purposes, we follow Jeffreys (1965). 


We first consider conventional Legendre polynomials, which can be defined 

with 
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2 i '. dz 
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For computational purposes, we can use 
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where we use {x} to be the greatest integer in x. These polynomials are orthogonal 
such that 


A-l 



f P„ (z) P fl , ,(z) r sin l mX r sin l m'X cos <j> d4> dX 
J , fm' ' i'm ,v ' |_cosJ LcosJ T 

sphere 


2tt (l + m) I 
e m (2* + l)(l - m)! 



= 0 


H. =M', m ^m', 
or both 


where 


_ ( 1 m = 0 
€ m (2 m + 0 


Each P im (z) can have a scale factor, called the normalization. We can choose this 
scale factor such that 
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We note that Jeffreys uses 
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We want the expression for p ^ m ( z ) e in another coordinate system defined by 
the conventional Euler angles (see Fig. Al). We can write 
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We can write this in a more compact form, if c|>' = 0, as 
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Reference : 

Jeffreys, B. S. , 1965. Transformation of tesseral harmonics under rotation. 
Geophys. Joum. Roy. Astron. Soc., vol. 10, pp. 141-145. 


A -4 



302-110 



Figure Al. Geometry of coordinate transformation. 
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